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1.  Introduction 


Vertical  profiles  of  meteorological  variables  produced  by  numerical  weather 
prediction  (NWP)  models  can  be  compared  to  data  from  World  Meteorological 
Organization  (WMO)  and  other  radiosondes.  Other  sources  of  comparison  data 
may  come  from,  for  example,  radar  profilers  and  lidar  for  wind,  microwave 
radiometers  for  temperature  and  moisture  parameters,  and  radio  acoustic  sounding 
systems  (RASSs)  for  virtual  temperature.  Nevertheless,  radiosondes  have  remained 
the  primary  source  of  comparison  data  above  the  near  surface  layer,  starting  roughly 
tens  of  meters  to  perhaps  100  m  above  the  surface. 

Meteorological  centers  compare  global  and  large-scale  regional  model  output  vs. 
radiosonde  data  more  or  less  continuously  and  make  intermodel  comparisons.  The 
National  Center  for  Environmental  Prediction  (NCEP)  has  a  readily  accessible 
website  (http://www.emc.ncep.noaa.gov/gmb/STATS  vsdb/)  and  the  European 
Center  for  Medium-Range  Weather  Forecasting  (ECMWF)  provides  related 
information  on  their  website  (http  ://www.  ecmwfi  int / en/  forecasts/ charts/ 
medium/monthly-wmo-scores-against-radiosondes).  Worldwide  comparisons  are 
available  for  deterministic  forecasts  at  the  ECMWF  (http://apps.ecmwf.int/ 
wmolcdnv/)  and  ensemble  forecasts  at  the  Japan  Meteorological  Agency  (JMA) 
(http://epsv.kishou.go.jp/EPSv/).  They  provide  comparisons  of  meteorological 
variables  at  the  standard  WMO  pressure  levels  (e.g.,  850,  700,  500,  400  hPa). 

Some  documentation  with  respect  to  these  comparisons  and  model  details  may  be 
found  via  the  respective  web  pages.  For  the  ECMWF,  one  can  go  to 
http://www.ecmwf.int/search/elibrary/.  For  the  NCEP,  the  following  page  and 
included  links  lead  to  various  documents  and  other  information  on  models  and 
related  datasets:  http://www.emc.ncep.noaa.gov/?doc=doc. 

Somewhat  similar  comparisons  of  radiosondes  and  smaller-scale  models  may  be 
found  in  the  published  literature.  Schroeder  et  al.  (2006)  investigated  most  of  the 
vertical  extent  of  the  Meteorological  Model  Fifth  Generation  (MM5)  and  presented 
results  for  standard  pressure  levels  from  850  to  100  hPa  (some  graphs  to  150  hPa) 
in  their  evaluation  of  their  automated  rapidly  relocatable  nowcasting  and  prediction 
system.  They  worked  with  data  for  8  days  in  April  2002,  12  days  in  the  winter  and 
summer  of  2003,  and  18  days  in  August  2001.  The  former  2  periods  were  from  the 
East  Coast  region  of  the  United  States  and  the  latter  from  the  Great  Plains  region 
of  the  United  States.  Cuevas  et  al.  (2011)  made  21  Weather  Research  and 
Forecasting  (WRF)  to  radiosonde  comparisons  in  Chile  during  24  October-4 
November  2011.  While  the  emphasis  in  the  report  concerned  precipitable  water 
vapor  (PWV)  forecasts  for  use  by  the  Southern  European  Observatory,  it  presented 
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vertical  profiles  of  the  mean,  bias,  and  root  mean  squared  error  (RMSE)  of  several 
other  meteorological  variables  such  as  temperature,  humidity,  and  wind  speed  that 
were  compared  as  a  function  of  pressure  level.  Cortes  and  Cure  (2011)  compared 
results  from  the  Global  Forecast  System  (GFS),  the  ECMWF  model,  and  WRF  at 
the  standard  pressure  levels  for  locations  in  northern  Chile.  Kilpelainen  et  al.  (2012) 
and  Dutsch  (2012)  each  evaluated  WRF  output  for  the  boundary  layer  over 
Svalbard  in  the  Arctic  in  tenns  of  height  above  ground  compared  to  tower  and 
tethered  balloon  (tethersonde)  data  and  radiosonde  data,  respectively.  The 
tethersondes  provided  data  from  the  surface  up  to  about  600,  800,  and  1250  m  for 
each  of  3  sites.  The  radiosonde  data  reached  as  high  as  2  km  above  the  surface. 
Additionally,  Behne  (2008)  addressed  wind  speed  as  derived  from  wind 
components  (vector  wind  speed/magnitude). 

Other  means  to  more  fully  evaluate  vertical  profiles  derived  from  model  output 
should  help  in  the  evaluation  of  model  output  as  compared  to  observations,  for 
example,  as  compared  to  data  from  radiosondes  or  for  wind  to  data  from  radar 
profilers  or  lidars.  Here  we  develop  methods  to  produce  vertical  profiles  of 
meteorological  variables  in  terms  of  height  and  pressure  levels,  and  generate 
integrated  mean  value  profiles  of  those  variables  for  layers  as  defined  by  the  user- 
input  values  of  boundary  heights  or  pressures.  The  output  profiles  are  entered  into 
spreadsheets  for  calculation  of  the  derived  variables  of  density  and  vector  wind 
speed  or  vector  wind  magnitude  where  those  latter  terms  refer  to  wind  speed 
derived  from  the  u  and  v  wind  components. 

The  level  and  mean  layer  values  from  WRF  (version  3.6.1)  are  compared  to  values 
from  co-incident  WMO  radiosonde  observation  (RAOB)  data  and  the  differences 
computed  for  the  several  variables.  Skamarock  et  al.  (2008)  describe  the  basics  of 
the  WRF  model,  though  there  have  been  some  changes  in  the  more  recent  versions. 
In  addition  statistical  values  are  produced  for  each  level  or  layer.  Though  the  set  of 
data  is  not  large  enough  to  be  definitive  with  respect  to  WRF  accuracy,  the  statistics 
appear  reasonable,  and  the  values  for  model  output  pressure  levels  or  equivalent 
heights  are  generally  in  line  with  previously  published  values  as  well  as  with  values 
for  pressure  levels  from  global  models  at  the  operational  sites  noted  above.  Here 
accuracy  is  defined  as  the  closeness  to  values  from  co-incident  radiosondes.  The 
methods  developed  here  may  be  used  for  a  variety  of  model  generated  and  observed 
vertical  profiles  with  only  minimal  changes,  primarily  to  the  function  that  reads  the 
input  data. 
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2.  Method 


The  generation  of  the  4  types  of  output  described  in  this  report  is  based  on  the 
algorithms  and  software  described  in  Cogan  (2015).  Here  we  describe  the  height- 
and  pressure-based  methods;  for  each  type,  there  is  one  method  for  values  at  user 
selected  levels  and  one  for  the  layers  defined  by  those  levels.  The  selected  heights 
or  pressure  levels  are  entered  via  a  user-generated  parameter  file.  It  contains  a  list 
of  the  desired  heights  or  pressures  starting  from  the  lowest  height  in  terms  of  above 
ground  level  (AGL)  or  the  highest  pressure,  respectively,  normally  values  for  the 
surface.  The  4  types  of  vertical  profiles  are  height  level,  height  layer,  pressure  level, 
and  pressure  layer.  These  methods  were  used  to  produce  RAOB-  and  WRF -based 
profile  comparisons  for  2  times  for  each  of  15  WMO  radiosonde  sites  in  the 
Southwest  United  States,  Mid-Atlantic  United  States,  Gennany,  and  South  Korea 
for  a  total  of  30  comparisons. 

2.1  RAOB  Soundings  and  WRF  Output  for  Profile  Generation 

The  RAOB  data  were  obtained  from  the  University  of  Wyoming’s  weather  website 
(http://www.weather.uwyo.edu/upperair/sounding.html).  That  site  contains  WMO 
soundings  in  several  formats  including  text  as  used  here.  Other  sites  are  available 
that  have  WMO  sounding  data  but  require  different  processing  in  the  input  function 
(e.g.,  http://www.esrl.noaa.gov/raobs).  The  sounding  data  are  transferred  to  text 
files  for  input  to  the  program  as  in  Cogan  (20 1 5).  As  in  the  earlier  program,  the  data 
are  entered  into  a  “standard”  array  for  further  processing.  One  input  function  is  used 
for  all  of  the  methods  for  radiosonde  sounding  data,  and  a  second  input  function  for 
all  of  those  for  model  output  profiles. 

WRF  v3.6.1  was  run  with  9/3/1  km  horizontal  grid  spacing  nested  domains.  The  30 
comparisons  for  this  report  used  data  from  the  3 -km  domain.  The  initial  and 
boundary  conditions  were  derived  from  GFS  0.5°  horizontal  grid  spacing  with  a 
3-h  time  interval.  Where  available,  GFS  snow  fields  were  replaced  with  1)  1  -km 
snow  fields  from  the  National  Weather  Service’s  National  Operational  Hydrologic 
Remote  Sensing  Center  (NOHRSC)  (http://www.nohrsc.noaa.gov/technology/) 
Snow  Data  Assimilation  System  (SNODAS),  or  if  1)  was  not  available  with 
2)  4-km  snow  cover  fields  from  the  National  Ice  Center’s  Interactive  Multi-sensor 
Snow  and  Ice  Mapping  System  (IMS)  (http://nsidc.org/data/docs/noaa/ 
g02156  ims  snow  ice  analysis/).  Consequently,  when  available,  NOHRSC  fields 
were  used  for  an  area  centered  on  the  United  States,  IMS  data  for  the  rest  of  the 
northern  hemisphere,  and  GFS  fields  for  the  southern  hemisphere. 
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A  sea  surface  temperature  product  with  higher  resolution  than  the  GFS  output  is 
produced  by  the  NCEP  Prediction,  Marine  Modeling  and  Analysis  Branch,  called 
the  Real  Time  Global  Sea  Surface  Temperature  (Gemmill  et  al.  2007),  which  has 
l/12th-degree  horizontal  grid  spacing  and  was  used  to  specify  sea  surface 
temperatures.  Observation  “nudging”  data  assimilation  ingested  observations 
during  a  3.5-h  preforecast  period.  Above-surface  observations  are  used  in 
assimilation  if  they  are  within  1.5  h  of  the  current  time,  while  the  time  window  for 
surface  observations  is  75%  as  long  (67.5  min).  It  nudged  the  model  toward 
observations  of  temperature,  moisture,  and  wind;  these  observations  when  available 
include  Meteorological  Assimilation  Data  Ingest  System  (MADIS)  data  from 
radiosondes,  Meteorological  Aerodrome  Report  (METAR),  Surface  Aviation 
Observation  (SAO,  mostly  from  Canada),  and  maritime  (e.g.,  ship)  observations. 
These  observations  are  only  entered  if  they  pass  the  quality  control  procedures  and 
are  within  the  respective  time  window.  The  Mellor-Yamada-Janjic  scheme  (MYJ) 
is  used  to  parameterize  the  atmospheric  boundary  layer.  As  in  Lee  et  al.  (2012)  and 
Reen  et  al.  (2014),  the  background  turbulent  kinetic  energy  is  decreased  to  better 
simulate  conditions  with  low  turbulent  kinetic  energy  and  the  atmospheric 
boundary  layer  depth  diagnosis  is  altered.  The  WRF  single-moment  5-class 
microphysics  parameterization  and  the  Kain-Fritsch  cumulus  parameterization 
(9-km  domain  only)  are  used.  For  radiation,  the  Rapid  Radiative  Transfer  Model  is 
used  for  longwave  and  the  Dudhia  scheme  for  shortwave.  The  Noah  land  surface 
model  is  used  to  represent  land  surface  processes.  The  pre-processing  software 
included  updates  by  Reen  (2015). 

WRF  was  run  for  21  model  hours  starting  3  h  after  either  00  UTC  or  12  UTC  with 
data  extracted  from  the  9-h  files  leading  to  output  for  12  UTC  or  00  UTC  the  next 
day,  respectively.  The  model  time  of  9  h  was  chosen  so  that  the  output  could  serve 
additional  investigations  beyond  the  scope  of  this  study.  The  profiles  of 
meteorological  variables  were  then  extracted  from  the  WRF  NetCDF  output  files 
via  a  National  Center  for  Atmospheric  Research  (NCAR)  Command  Language 
(NCL)  script  and  placed  into  text  files  that  included  a  header  with  infonnation  on 
the  location,  time,  model  grid  resolution,  and  the  method  (if  any)  employed  to 
interpolate  between  model  grid  points  to  the  selected  location  (Reen  2015).  These 
profiles  contain  data  lines  for  heights  above  mean  sea  level  (MSL)  of  pressure 
levels.  The  user  has  the  option  to  use  a  “sounding”  from  the  nearest  grid  point  (i.e., 
no  interpolation),  or  either  bilinear  interpolation  or  inverse  distance  weight 
interpolation  (via  built-in  NCL  capabilities). 
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2.2  Height-Based  Profiles 


The  software  for  height-based  output  required  minimal  modification  to  the  program 
described  in  Cogan  (2015)  including  changes  to  produce  height  output  at  specific 
heights  as  well  as  layer  output.  This  revised  version  also  is  written  in  the  “C” 
language  and  that  language’s  terminology  is  used  herein,  such  as  the  term  function 
for  a  sub-program  or  routine.  Appendix  A  contains  details  on  the  changes  to  the 
earlier  software,  but  does  not  repeat  the  entire  program. 

The  changes  to  the  earlier  software  mostly  involved  the  addition  of  output  at  the 
user-defined  heights.  Those  values  were  computed  previously  as  part  of  the 
calculation  of  the  mean  layer  values,  but  were  not  included  in  the  output.  The 
present  version  has  an  additional  function  that  produces  height  level  output.  As 
before,  the  user  defines  the  output  height  profile  in  a  parameter  file  and  that,  in  turn, 
defines  the  upper  and  lower  boundaries  of  the  computed  integrated  mean  layer 
values  (a.k.a.  weighted  mean  values).  For  both  the  level  and  layer  output,  the  initial 
data  line  contains  the  input  profile’s  surface  values,  but  converted  to  the  output 
fonnat.  Subsequent  data  lines  denote  either  level  or  layer  output,  respectively.  For 
the  level  output,  the  values  displayed  are  those  computed  for  the  listed  heights.  The 
layer  “height”  displayed  is  the  midpoint  height  of  the  layer,  but  a  given  data  value 
may  not  represent  the  exact  value  at  the  midpoint  height  since  it  is  the  weighted 
mean  of  that  layer.  The  exception  is  pressure,  which  is  calculated  for  the  midpoint 
using  the  hypsometric  formula  starting  from  the  closest  data  level  immediately 
below  or  the  lower  boundary  value  if  there  is  no  data  level  between  the  lower 
boundary  height  and  the  midpoint. 

The  input  list  of  heights  in  the  parameter  file  may  be  modified  by  the  user  as  long 
as  there  are  2  levels,  a  “surface”  (at  the  surface  or  defined  lowest  level)  and  top 
level.  Heights  are  listed  as  AGL.  With  minimal  modification  the  program  will  list 
heights  as  above  MSL  in  the  output.  The  user-input  levels  may  define  a  vertical 
spacing  or  interval  between  heights  or  layer  thickness  from  a  fraction  of  meter  to 
the  entire  vertical  extent  that  can  exceed  several  kilometers.  However,  the  user 
should  be  aware  of  the  vertical  spacing  between  input  data  levels  of  the  RAOB  or 
extracted  WRF  profiles  before  using  a  very  small  or  large  interval. 

2.3  Pressure-Based  Profiles 

The  methods  presented  in  this  section  also  are  based  on  those  presented  in  Cogan 
(2015)  but  required  more  extensive  modification  than  for  height -based  output. 
Appendix  A  contains  details  on  these  changes  but,  as  with  the  height-based 
versions,  does  not  repeat  the  entire  program.  Earlier  and  ongoing  efforts  as 
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referenced  in  the  introduction  provide  comparisons  for  standard  pressure  levels  or 
model  output  pressure  levels  for  specified  periods  during  field  tests.  The  methods 
derived  herein  are  not  limited  to  either  standard  or  model  output  levels.  As  with  the 
height-based  methods,  the  user  defines  the  vertical  levels,  in  this  case,  pressure 
instead  of  height. 

The  software  produces  profiles  of  the  several  variables  for  pressure  levels  or  layers 
specified  via  the  user’s  parameter  file.  The  “surface”  values  for  either  the  level  or 
layer  output  repeat  the  input  surface  values,  and  that  input  is  used  to  compute  the 
derived  surface  variables  (e.g.,  density),  with  appropriate  changes  in  units.  The  user 
may  specify  pressure  levels  that  produce  intervals  that  range  from  a  fraction  of  a 
hectopascal  to  the  entire  vertical  extent,  but  the  user  has  to  include  at  least  a  lowest 
pressure  level  (highest  pressure)  and  the  top  pressure  level  (lowest  pressure). 
However,  as  a  general  rule  the  pressure  level  intervals  should  not  be  too  large  (e.g., 
not  greater  than  about  50  hPa)  and  consider  the  altitude  since  a  change  of,  say,  25 
hPa,  is  not  too  large  near  the  surface,  but  may  be  near  the  50  hPa  level.  Also,  the 
input  pressure  intervals  should  be  considered  since  a  very  fine  output  interval  (e.g., 
1  hPa  or  less)  may  not  yield  much  additional  infonnation  if  the  input  interval  is 
much  larger  (e.g.,  25  hPa). 

Modifications  to  the  software  to  enable  pressure-based  output  included  changes  to 
the  output  functions  and  modification  of  the  function  that  produces  the  level  and 
layer  values.  The  standard  version  of  the  program  linearly  interpolates  all  variables 
(e.g.,  height,  temperature)  in  terms  of  the  natural  log  of  pressure,  ln(P).  An 
alternative  version  of  the  program  computes  height  for  a  user-defined  level  using 
the  hypsometric  formula  from  the  input  pressure  level  immediately  below  (next 
highest  input  pressure);  this  alternative  version  retains  linear  in  ln(P)  interpolation 
for  all  other  variables.  If  no  data  level  lies  between  the  given  data  level  and  the 
output  data  level  immediately  below,  then  the  program  uses  the  value  of  pressure 
at  that  previous  output  level.  The  same  procedure  applies  to  computation  of  the 
layer  value  of  height  except  that  the  hypsometric  fonnula  is  applied  to  the  midpoint 
of  the  data  layer  (Fig.  1).  Figure  2  illustrates  the  computation  of  height  using  linear 
interpolation  in  ln(P). 
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Upper  boundary  level 


o 


Midpoint 


Lower  boundary  level 


Fig.  1  Schematic  chart  based  on  a  similar  diagram  in  Cogan  (2015)  illustrating  the  ln(P) 
level  from  which  the  height,  Z,  of  the  midpoint  is  computed  using  the  hypsometric  equation. 
“A”  represents  computation  starting  at  the  data  point  immediately  below  the  midpoint  (solid 
circle),  and  “B”  indicates  computation  starting  at  the  lower  boundary  level  (solid  square)  of  a 
given  layer  (only  undertaken  if  the  point  used  in  “A”  is  not  present). 


Upper  boundary  level 


Midpoint 


Lower  boundary  level 


Fig.  2  Schematic  chart  illustrating  the  linear  interpolation  of  height  in  terms  of  ln(P)  for  2 
cases  where  the  solid  triangles  represent  the  interpolated  values.  The  first  (A)  is  an  example 
of  interpolation  between  2  data  points  that  lie  within  the  layer  (solid  circles).  The  second  (B) 
is  an  example  of  interpolation  between  boundary  levels  (solid  squares)  where  there  are  no  data 
levels  within  the  layer.  Intermediate  cases  can  occur  where  a  data  level  is  within  the  layer 
above  or  below  the  midpoint. 

Several  comparisons  of  the  2  methods  for  height  showed  that  the  output  differences 
in  height  of  pressure  levels  were  small  (absolute  difference  of  0  to  3  m  for  most 
RAOB  levels  to  the  nearest  meter,  and  0  to  2  m  for  WRF  levels).  The  differences 
in  the  layer  values  were  somewhat  larger,  especially  for  RAOB  pressures  less  than 
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40  hPa.  Standard  WRF  procedure  allows  interpolation  in  either  linear  or  linear  in 
natural  log  of  pressure  (via  option  interp  type  as  noted  in  http:// 
www2.mmm.ucar.edu/wrf/users/docs/user  guide  V3/users  guide  chap5.htm  ), 

The  primary  temperature  variable  in  WRF  is  potential  temperature  and  WRF 
interpolates  it  in  terms  of  ln(P)  even  if  the  user  selects  linear  in  P  for  the  other 
variables. 

The  procedure  for  computation  of  pressure  layer  values  closely  resembles  that  for 
computation  of  height  layer  values.  The  main  difference  is  that  the  vertical 
coordinate  is  ln(P)  instead  of  height  (Z).  Figure  3  illustrates  the  method  for  either 
Z  or  ln(P)  for  any  of  the  variables,  except  for  P  in  the  height-based  computations 
where  the  hypsometric  fonnula  is  used  to  compute  P  as  well  as  Z  in  the  alternate 
pressure  method  that  also  uses  the  hypsometric  fonnula. 


Fig.  3  Schematic  based  on  a  similar  diagram  in  Cogan  (2015)  illustrating  the  computation 
of  an  integrated  or  weighted  mean  of  some  variable  X  for  a  layer  defined  by  height  (Z),  upward 
arrow,  or  the  natural  log  of  pressure,  ln(P),  downward  arrow.  The  blue  lines  represent 
sublayers  with  sublayer  means  indicated  by  red  triangles.  Circles  indicate  the  observations  or 
WRF  output  values  and  blue  squares  the  values  at  the  upper  and  lower  layer  boundary  levels. 


3.  Comparisons 

For  each  method,  the  data  for  each  pair  of  RAOB  and  WRF  computed  level  and 
each  pair  of  layer  profile  text  files  were  entered  into  separate  spreadsheets  in  the 
form  of  tables  by  variable  and  level  or  layer.  Table  1  provides  an  example  of 
computed  height  level  values  from  the  00  UTC  RAOB  from  Wallops  Island, 
Virginia,  on  27  October  2015,  which  has  the  same  format  as  the  co-incident  data 
from  WRF.  The  heights  are  in  meters  above  ground  level  (m  AGL).  Only  the  first 
12  of  60  levels  are  shown. 
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Table  1  Sample  of  height  level  “sounding”  derived  from  RAOB  data  at  00  UTC  on  27 
October  2015  at  Wallops  Island,  Virginia.  Height  is  in  m  AGL,  wind  direction  (W  dir)  in 
degrees,  wind  speed  (W  spd)  in  knots,  virtual  temperature  (Tv)  and  sensible  temperature  (T) 
in  K,  and  pressure  (P)  in  hPa. 


Level 

Height 

W  dir 

W  spd 

Tv 

P 

T 

0 

0 

350 

6.0 

285.7 

1029.0 

284.6 

1 

50 

6 

7.2 

285.5 

1022.9 

284.5 

2 

100 

16 

8.7 

285.4 

1016.8 

284.4 

3 

200 

38 

11.8 

284.9 

1004.3 

284.0 

4 

300 

50 

16.0 

284.1 

992.1 

283.2 

5 

400 

54 

16.3 

283.2 

980.2 

282.3 

6 

500 

57 

16.6 

282.2 

968.4 

281.4 

7 

600 

60 

17.0 

281.3 

956.6 

280.4 

8 

700 

64 

16.6 

280.4 

945.0 

279.5 

9 

800 

67 

16.2 

279.5 

933.5 

278.6 

10 

900 

70 

16.0 

278.6 

922.4 

277.8 

11 

1000 

63 

16.0 

278.1 

911.1 

277.4 

Table  2  presents  a  similar  example,  but  for  computed  pressure  level  values  from 
the  00  UTC  RAOB,  also  from  Wallops  Island  on  27  October  2015.  The  co-incident 
WRF  “sounding”  has  the  same  format  and  pressure  levels.  The  pressure  values 
shown  are  those  listed  in  the  user-provided  parameter  file  and  almost  always  are 
above  the  surface.  An  exception  could  occur  if  both  the  WRF  and  RAOB  profiles 
had  the  same  surface  pressure  and  that  was  the  same  as  a  value  in  the  parameter 
file. 


Table  2  Sample  of  pressure  level  “sounding”  derived  from  RAOB  data.  Pressure  (P)  is  in 
hPa,  height  in  m  AGL,  wind  direction  (W  dir)  in  degrees,  wind  speed  (W  spd)  in  knots,  and 
virtual  temperature  (Tv)  and  sensible  temperature  (T)  in  K.  Data  for  Wallops  Island,  Virginia, 
on  27  October  2015  at  00  UTC. 


Level 

P 

Height 

W  dir 

W  spd 

Tv 

T 

0 

1020 

72 

11 

7.8 

285.5 

284.5 

1 

1010 

152 

27 

10.4 

285.2 

284.3 

2 

1000 

234 

45 

13.0 

284.7 

283.8 

3 

990 

317 

51 

16.1 

283.9 

283.0 

4 

975 

443 

55 

16.4 

282.8 

281.9 

5 

950 

657 

62 

16.7 

280.8 

279.9 

6 

925 

877 

70 

16.0 

278.8 

278.0 

7 

900 

1100 

57 

16.0 

280.5 

280.4 

8 

875 

1332 

52 

14.5 

281.4 

281.4 

9 

850 

1571 

55 

11.0 

280.6 

280.6 

10 

825 

1815 

50 

11.0 

279.8 

279.8 

11 

800 

2067 

1 

5.1 

279.0 

279.0 

The  differences  between  the  WRF  and  RAOB  values  are  listed  in  a  table  on  the 
same  spreadsheet  by  each  variable  and  level  or  layer.  In  addition  to  the 
aforementioned  variables,  spreadsheet  functions  compute  density  (Den)  difference, 
in  gm  3  and  vector  wind  magnitude  (V-W)  difference,  in  knots.  The  difference 


Approved  for  public  release;  distribution  unlimited. 


9 


values  for  each  variable  is  written  to  a  summary  spreadsheet  for  each  of  the  levels 
or  layers.  The  summary  sheet  is  set  up  so  that  standard  statistics,  i.e.,  mean  value 
(M),  mean  absolute  error  (MAE),  standard  deviation  (SD),  and  RMSE,  are 
computed  for  each  level  or  layer  for  each  variable.  For  this  set  of  comparisons 
“error”  really  refers  to  difference,  and  from  here  on,  this  report  uses  the  terms  mean 
difference  (MD),  mean  absolute  difference  (MAD),  SD  of  the  differences,  and  root 
mean  square  difference  (RMSD).  Table  3  presents  a  sample  of  the  height-based 
difference  output  for  the  same  set  of  profiles,  that  is,  for  Wallops  Island  on  27 
October  2015  at  00  UTC.  Some  intennediate  values  have  been  omitted  such  as  the 
outcomes  of  modifying  the  wind  direction  by  adding  360°  to  one  of  the  profiles  if 
the  absolute  value  of  the  initial  direction  difference  would  exceed  180°  (e.g.,  if 
WRF  value  -  RAOB  value  =  1 10  —  350|  =  340°,  add  360  to  the  WRF  value:  10+360 
-  350  =  20°  difference).  Only  the  first  12  of  60  levels  are  shown. 

Table  3  Sample  of  differences  in  listed  variables  derived  from  WRF  and  RAOB  data 
(WRF  -  RAOB)  for  listed  height  (m  AGL)  levels.  Difference  is  in  wind  direction  (W  dir)  is  in 
degrees,  wind  speed  (W  spd)  in  knots,  virtual  temperature  (Tv)  and  sensible  temperature  (T) 
in  K,  pressure  (P)  in  hPa,  density  (Den)  in  gm-3,  and  vector  wind  magnitude  (V-W)  in  knots. 
Data  for  Wallops  Island,  Virginia,  on  27  October  2015  at  00  UTC. 


Level 

Height 

W  dir 

W  spd 

Tv 

T 

P 

Den 

V-W 

0 

0 

55 

6.3 

1.1 

1.1 

0.7 

-3.96 

10.13 

1 

50 

41 

8.0 

1.6 

1.4 

0.5 

-6.35 

10.85 

2 

100 

33 

8.8 

1.4 

1.2 

0.6 

-5.33 

11.25 

3 

200 

13 

7.3 

1.0 

0.8 

1.0 

-3.08 

8.05 

4 

300 

2 

3.7 

0.8 

0.6 

1.2 

-1.95 

3.75 

5 

400 

-1 

3.6 

0.7 

0.5 

1.2 

-1.50 

3.61 

6 

500 

-4 

3.4 

0.7 

0.4 

1.3 

-1.36 

3.63 

7 

600 

-6 

3.1 

0.6 

0.5 

1.4 

-0.79 

3.65 

8 

700 

-10 

3.2 

0.6 

0.5 

1.5 

-0.65 

4.50 

9 

800 

-9 

2.1 

1.0 

1.3 

1.5 

-2.29 

3.42 

10 

900 

-12 

-0.1 

2.3 

2.8 

1.3 

-7.83 

3.34 

11 

1000 

-5 

-2.0 

3.3 

3.9 

1.4 

-11.65 

2.39 

A  similar  table  of  difference  values  was  computed  for  the  pressure  level-based  data, 
but  with  a  column  for  height  differences  vs.  pressure  differences.  The  same  types 
of  tables  were  computed  for  the  layer  output.  However,  for  the  height-based  data, 
the  listed  heights  were  the  layer  midpoints.  For  the  pressure -based  data,  the  listed 
pressures  were  the  averages  of  the  upper  and  lower  boundary  pressures.  Table  4 
presents  difference  data  from  Kwangju,  South  Korea,  for  15  September  2015  at 
12  UTC  for  the  surface  and  first  6  of  59  height  layers.  Level  0  contains  the  surface 
data  with  format  changes  and  variables  derived  from  those  data  (e.g.,  density).  The 
layer  0  data  line  has  the  same  values  as  the  height  level  0  line.  Layers  starting  at  1 
contain  layer  values. 
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Table  4  Difference  data  as  in  Table  3  (WRF  -  RAOB),  but  for  the  surface  and  first  6  height 
layers.  The  variables  have  the  same  units  as  in  Table  3.  The  listed  heights  except  for  the  surface 
(laver  0)  are  the  midpoints  of  the  respective  lavers.  Data  for  Kwangju,  South  Korea,  for  15 
September  2015  at  12  UTC. 


Layer 

Fleight 

W  dir 

Wspd 

Tv 

T 

P 

Den 

V-W 

0 

0 

52 

4.0 

-3.9 

-3.9 

1.2 

17.27 

4.00 

1 

25 

-20 

4.3 

-2.6 

-2.5 

1.1 

11.80 

4.35 

2 

75 

-2 

5.7 

-1.3 

-1.3 

1.0 

6.38 

5.70 

3 

150 

3 

5.7 

-1.3 

-1.3 

0.6 

5.88 

5.71 

4 

250 

-7 

4.9 

-1.3 

-1.5 

0.6 

5.85 

5.03 

5 

350 

-13 

3.3 

-1.4 

-1.6 

0.6 

6.22 

4.19 

6 

450 

-16 

1.8 

-1.5 

-1.7 

0.5 

6.47 

4.23 

Table  5  presents  difference  data  from  Kwangju,  South  Korea,  for  15  September 
2015  at  12  UTC  for  the  first  6  of  62  pressure  layers.  For  pressure -based  data,  if  the 
highest  pressures  for  the  WRF-  and  RAOB-based  profiles  are  above  the  specified 
pressure  of  the  lower  boundary  of  the  initial  user-defined  layer  then  the  output  will 
include  that  layer  (layer  0).  Otherwise,  the  first  user  layer  to  be  output  is  the  one 
immediately  above  (i.e.,  layer  1)  and  so  on  for  subsequent  layers.  In  the  example 
shown,  the  highest  pressure  for  the  WRF  as  well  as  RAOB  profiles  were  less  than 
the  1020  hPa  specified  for  the  lower  boundary  of  layer  0  and  therefore  no  1015  hPa 
layer  difference  values  were  computed. 

Table  5  Difference  data  (WRF  -  RAOB)  as  in  Table  4,  but  for  the  first  6  pressure  layers. 
The  variables  have  the  same  units  as  in  Table  4.  The  listed  pressures  are  the  “midpoints” 
(average  of  boundary  pressures)  of  the  respective  lavers.  Data  for  Kwangju,  South  Korea,  for 
15  September  2015  at  12  UTC. 


Layer 

Press 

W  dir 

Wspd 

Tv 

T 

P 

Den 

V-W 

1 

1005.0 

3 

5.9 

-1.3 

-1.3 

7.0 

5.20 

5.91 

2 

995.0 

0 

5.7 

-1.4 

-1.4 

6.0 

5.56 

5.70 

3 

982.5 

-9 

4.4 

-1.4 

-1.6 
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Several  standard  statistical  measures  were  computed  for  each  of  the  4  methods  for 
each  variable  (direct  or  derived)  for  all  of  the  individual  levels  or  layers  covering 
all  30  radiosondes  used  for  comparison  against  WRF  model  output.  As  noted  above 
the  measures  were  MD,  MAD,  SD,  and  RMSD.  Tables  of  the  statistical  values  were 
generated  and  entered  into  tables  on  a  spreadsheet,  which,  in  turn,  were  used  to 
produce  graphical  representations.  Figures  4-7  present  those  4  statistics  for  height 
levels  and  layers  and  pressure  levels  and  layers  for  temperature  and  vector  wind 
magnitude.  The  vector  wind  magnitude  is  always  positive,  hence  the  MD  and  MAD 
values  are  the  same.  On  the  scale  of  the  charts,  SD  and  RMSD  frequently  overlap. 
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The  complete  set  of  charts  for  the  height  and  pressure  levels  and  layers  measures 
for  all  variables  are  presented  in  Appendix  B. 
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Fig.  4  Statistical  values  (MD,  MAD,  SD,  and  RMSD)  for  vector  wind  magnitude  and 
temperature  for  the  user-defined  height  levels  along  with  the  number  of  samples  for  each  level 
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Fig.  5  Statistical  values  (MD,  MAD,  SD,  and  RMSD)  for  vector  wind  magnitude  and 
temperature  for  height  layers  defined  by  the  user  input  height  levels  along  with  the  number 
of  samples  for  each  layer 
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Vector  Wind  Speed  Difference 


Fig.  6  Statistical  values  (MD,  MAD,  SD,  and  RMSD)  for  vector  wind  magnitude  and 
temperature  for  pressure  levels  defined  by  the  user  along  with  the  number  of  samples  for  each 
level 
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Vector  Wind  Magnitude  Difference 
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Fig.  7  Statistical  values  (MD,  MAD,  SD,  and  RMSD)  for  vector  wind  magnitude  and 
temperature  for  pressure  layers  defined  by  the  user-defined  levels  along  with  the  number  of 
samples  for  each  layer 


Earlier  published  and  ongoing  comparisons  at  operational  centers  as  noted  in  the 
references  and  websites  in  Section  1  provide  values  of  some  meteorological 
variables  between  model  output  and  radiosonde  soundings  for  several  standard 
pressure  levels  and  some  limited  comparisons  for  model  output  pressure  levels  and 
equivalent  heights.  Converting  the  wind  units  here  from  knots  to  meter  per  second 
(m  s-1)  shows  that  the  results  here  are  roughly  comparable  with  respect  to  wind. 
That  also  holds  for  the  other  directly  measured  meteorological  variables.  The 
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computed  values  of  height  for  the  compared  standard  levels  (850,  500,  250  hPa)  on 
the  ECMWF  website  for  deterministic  output  at  24  h  of  model  time  for  the  northern 
hemisphere  from  several  national  models  had  noticeably  lower  values  for  850  and 
250  hPa  and  somewhat  lower  at  500  hPa. 

The  initialization  model  for  this  study,  GFS,  has  an  accuracy  close  to  or  about  equal 
to  that  for  the  other  models  based  on  the  data  presented  on  the  ECMWF  site 
(excluding  an  outlier  that  showed  greater  differences).  However,  the  locations 
examined  for  this  report  included  some  in  or  very  near  complex  terrain.  The  profile 
data  for  the  2  sites  with  the  largest  height  differences,  Kwangju,  Korea,  and 
Blacksburg,  Virginia,  had  the  WRF  surface  height  significantly  different  from  that  of 
the  radiosonde  launch  site,  about  35-38  m  lower.  The  respective  pressure  differences 
were  about  6  and  2  hPa  higher.  With  those  removed  the  computed  heights  were 
closer,  with  the  500  hPa  value  slightly  better  than  the  ECMWF  mean  value,  the  value 
for  850  hPa  was  slightly  higher  than  the  mean  of  the  models,  and  the  250  hPa  value 
was  within  the  spread  of  the  several  models  displayed.  The  discrepancies  with  respect 
to  surface  height  and  pressure  suggests  that  finer-scale  versions  of  WRF  (e.g.,  1  or 
0.5  km)  could  help  if  the  terrain  data  base  also  has  a  finer  grid  resolution  than  as  used 
in  this  study.  However,  that  remains  to  be  detennined. 

4.  Summary  and  Conclusion 

Two  methods  were  developed  to  provide  vertical  profiles  of  meteorological 
variables  and  some  derived  quantities  from  radiosonde  and  model  output  at  user- 
specified  height  or  pressure  levels.  Two  other  methods  provide  profiles  of 
integrated  or  weighted  mean  values  for  layers  based  on  those  user-specified  levels. 
Ongoing  comparisons  at  data  centers  have  concentrated  on  comparing  models  and 
some  of  those  models  with  radiosonde  data  at  standard  pressure  levels  plus  surface 
values.  Fimited  comparisons  elsewhere  for  soundings  derived  from  model  output 
for  data  levels  at  model  computed  pressure  levels,  or  equivalent  heights  (often 
within  or  not  far  above  the  boundary  layer),  have  occurred  for  relatively  short 
experimental  periods  at  specific  sites. 

The  emphasis  at  operational  data  centers  tends  toward  wind  speed,  temperature,  and 
geopotential  height,  though  other  variables  are  addressed  such  as  relative  humidity. 
The  methods  developed  here  can  produce  output  that  is  not  limited  to  standard 
pressure  levels  or  model  output  pressure  levels  (or  equivalent  heights).  Furthermore, 
these  methods  may  be  applied  to  applications  other  than  model  comparisons  with 
radiosonde  data  such  as  inter-comparison  of  observation  systems  such  as  radar  wind 
profilers  with  lidars,  comparisons  of  different  model  configurations,  and 
determination  of  the  best  types  of  vertical  profiles  for  specific  applications. 
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The  use  of  pressure-level  data  over  varying  terrain  may  cause  some  confusion  since 
for  a  given  level  some  locations  may  be  close  to  the  surface  while  others  are  not. 
For  example,  a  common  surface  pressure  at  Flagstaff,  Arizona  (2192  m  above  mean 
sea  level),  may  range  from  780  to  790  hPa.  Comparing  data  from,  say  Phoenix, 
Arizona  (384  m  above  mean  sea  level),  with  those  from  Flagstaff  for  the  850  hPa 
pressure  level  does  not  have  much  value.  However,  looking  at  data  for  a  comparison 
at  a  height  of,  say,  500  m  AGL  could  have  some  use  for  a  boundary  layer 
investigation.  This  “discrepancy”  between  common  practice  and  potential  needs 
suggests  a  hybrid  vertical  profde  to  compare  model  output  with  observations  above 
the  surface.  This  procedure  would  use  a  vertical  profde  based  on  height  AGL  for 
the  lower  part  of  a  hybrid  profile,  for  example,  from  the  surface  to  perhaps  the 
lowest  hundreds  of  meters  to  a  few  km  above  the  surface,  and  a  vertical  profile 
based  on  pressure  level  above.  The  same  could  apply  to  vertical  layer  profiles. 

The  use  of  vertical  profiles  of  atmospheric  layers  provides  a  means  to  sample  the 
atmosphere  in  a  semi-continuous  fashion  vs.  a  series  of  point  values  upward  from 
the  surface.  The  layers  can  be  very  thin  for  close  examination  and  allowed  to  be 
thicker  for  regions  of  less  detailed  study.  One  could  select  layers  only  10,  5,  or  even 
1  m  thick  near  the  top  of  the  boundary  layer  and  perhaps  500  m  thick  in  the  upper 
troposphere  before  approaching  the  tropopause  where  perhaps  50-m  layers  could 
be  useful.  If  pressure-level  layers  were  preferred  then  perhaps  1  to  50  hPa  layers 
could  be  appropriate,  depending  on  the  atmospheric  situation.  The  user  can  adjust 
the  vertical  resolution  of  the  output  profiles  via  modification  of  a  single  parameter 
file,  no  recompilation  or  other  adjustment  is  required. 

These  methods  were  applied  to  a  set  of  co-located  WRF  generated  profiles  and 
WMO  radiosonde  soundings.  They  provided  a  sample  of  the  use  of  the  4  methods 
for  comparing  profiles  from  WRF  model  output  with  those  from  RAOBs  for  height 
levels  and  layers  and  pressure  levels  and  layers.  Several  common  statistical 
measures  were  computed  from  the  resultant  difference  data.  Though  somewhat 
varied  in  elevation,  climate  region,  etc.,  the  30  cases  of  this  study  can  only  suggest 
the  potential  accuracy  of  WRF  output  relative  to  RAOBs  since  they  are  too  few  for 
a  definitive  conclusion  in  that  regard.  Here  accuracy  refers  to  closeness  to  the 
radiosonde  observations,  and  furthermore,  a  RAOB  is  only  a  good  estimate  of 
atmospheric  conditions  due  to  spatial  and  temporal  drift  and  instrument  errors. 
Nevertheless,  the  statistics  for  standard  pressure  levels  may  be  compared  with  those 
from  published  site-specific  results  or  displayed  on  public  websites  run  by 
meteorological  centers,  as  well  as  with  published  site-specific  results  for  model 
output  pressure  levels  or  equivalent  heights.  Overall  the  results  from  this  study  are 
in  line  with  the  values  presented  in  those  references  and  websites. 
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Appendix  A.  Program  Modification  Descriptions 
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This  appendix  contains  short  summaries  in  text  and  pseudo  code,  if  applicable,  of 
the  changed  sections  of  the  relevant  functions  in  the  programs  for  generation  of 
vertical  profdes  of  meteorological  variables  for  user-defined  heights,  pressure 
levels,  height  layers,  and  pressure  layers.  The  changes  to  the  height  based  methods 
were  relatively  small  since  the  generation  of  the  output  profiles  was  essentially  the 
same  as  in  Cogan. 1  There  were  output  fonnat  changes  and  other  very  minor 
modifications  aside  from  the  modification  to  the  parameter  file  that  contained  user- 
defined  heights.  Changes  to  the  pressure-based  methods  were  more  extensive,  but 
not  too  large.  The  output  function  required  some  modification  and  the  parameter 
file  contains  user-defined  pressure  levels  versus  heights.  A  more  extensive 
modification  concerned  the  computation  of  pressure  level  and  layer  values  of  the 
variables  as  well  as  the  calculation  of  height  level  and  layer  values  versus  pressure 
values. 

A-l  Height-Based  Methods 

The  main  changes  were  modifications  to  the  output  function  for  layer  data  and  a 
new  function  for  output  of  level  data,  plus  changes  in  fonnat.  For  example 
temperature  and  virtual  temperature  were  changed  from  K*  10  to  K.  Pressure  is  now 
labeled  in  hectopascal  (hPa)  instead  of  millibar  (mb)  to  comply  with  current 
common  usage  in  the  meteorological  community.  The  main  computation  function, 
“msgvalues”,  was  essentially  unchanged  from  that  described  in  Cogan.1  Since  the 
earlier  program  calculated  values  of  the  variables  for  height  levels  as  part  of  the 
computation  of  layer  values  for  the  several  types  of  MET  messages,  the  main 
change  was  to  retain  the  arrays  of  level  values  for  later  output.  A  separate  output 
function  was  created  for  the  level  output.  This  level  output  function  looks  almost 
the  same  as  the  layer  output  function  in  Cogan,1  but  the  heights  and  the  respective 
variables  are  for  the  levels  themselves  vs.  the  upper  boundary  heights  and  integrated 
mean  values  of  the  layers.  The  layer  output  function  now  writes  the  midpoint 
heights  of  the  layers  instead  of  the  heights  of  the  layer  upper  boundaries,  otherwise 
it’s  much  the  same  as  before.  The  changes  for  both  the  level  and  layer  output 
functions  were  very  small  and  consequently  are  not  shown  here.  Computation  of 
vector  wind  magnitude  and  density  were  performed  on  the  spreadsheets  using 
common  fonnulas  with  the  “direct”  profile  data  as  input. 


1  Cogan  J.  A  generalized  method  for  vertical  profiles  of  mean  layer  values  of  meteorological  variables. 
Adelphi  (MD):  Army  Research  Laboratory  (US);  2015.  Report  No.:  ARL-TR-7434. 
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A-2  Pressure-Based  Methods 


The  methods  for  pressure  (P)  levels  and  layers  required  greater  changes  to  the 
respective  programs  though  those  changes  were  not  large.  The  parameter  file  for 
user  input  contains  pressure  levels  versus  height  levels.  The  output  functions  and 
files  are  nearly  the  same  as  those  for  the  height-based  output.  Instead  of  output  of 
height  levels  or  height  layer  midpoints,  the  output  files  have  pressure  levels  or 
pressure  layer  midpoints.  Height  or  layer  value  heights  are  calculated  as  are  the 
other  meteorological  variables.  As  noted  in  Section  2  the  interpolation  is  in  terms 
of  the  natural  logarithm  of  P  (ln[P]),  not  P.  Consequently,  the  layer  values  of  height 
shown  are  not  strictly  the  values  at  the  P  midpoints,  but  rather  the  values  at  the  ln(P) 
midpoints.  A  limited  comparison  using  data  from  a  few  soundings  from  Dulles 
Airport,  Virginia,  and  Idar-Oberstein,  Gennany,  of  P  and  ln(P)  interpolation 
suggested  the  absolute  values  of  the  differences  in  computed  values  of  height  can 
be  small  (<1  m  up  to  4  m)  for  most  layers,  but  for  some  it  can  rise  to  over  10  m  and 
for  a  very  few  layers  in  the  radiosonde-based  profiles  to  over  20  m.  At  the  highest 
levels  in  the  radiosonde-based  profiles  (P  <  20  hPa)  where  temperature  increased 
with  height,  even  larger  differences  appeared  (maximum  magnitude  of  86  m  [P  vs. 
ln(P)  interpolated  heights]  for  the  10-15  hPa  layer  for  Dulles  Airport).  The  WRF- 
based  profiles  did  not  have  a  pressure  less  than  50  hPa,  or  in  other  words,  it  did  not 
go  above  about  20  km.  In  the  limited  set  of  data  examined  the  maximum  absolute 
values  of  level  or  layer  height  difference  was  7  m. 

A-2.1  Msgvaluesprs  Function 

Here  we  present  the  main  computation  function  (msgvaluesprs)  where  most  of  the 
substantive  changes  occurred.  This  function,  in  turn,  uses  2  generalized  functions, 
levelprs  and  layerprs,  that  perfonn  linear  interpolation  and  generation  of 
integrated/weighted  mean  layer  values  respectively  based  on  pressure  (P)  or  some 
function  of  pressure  such  as  ln(P)  as  was  done  here,  where  p  =  ln(P).  The 
description  of  msgvaluesprs  in  pseudo  code  and  text  follows.  The  pseudo  code  uses 
some  ‘C’  computer  language  syntax  such  as  /*  and  */  to  enclose  comments,  and  a 
semicolon  (;)  to  indicate  the  end  of  a  line  of  code.  The  natural  logarithm  in  ‘C’  is 
written  as  log  vs.  In.  The  parameter  ERROR  (=  -999)  indicates  missing  or  out  of 
bounds  data  and  is  used  to  initialize  most  arrays. 

Pass  in  the  input  data  structure  and  the  output  data  structure  for  the  output  layer  values  and  the  one 
for  the  computed  level  values. 

/*  Note  that  the  structures  are  structures  of  arrays  that  hold  header  information  as  well  as  the  data. 
Each  level  or  layer  is  an  element  of  an  array  that  contains  the  values  for  that  level  or  layer  as 
structure  elements.  Here  p  refers  to  ln(P).*/ 

{ 

/*  message  values  =  layer  values,  but  may  have  different  indices,  mlevel  =  level  values  */ 

Set  the  indices  and  the  variable  definitions. 


Approved  for  public  release;  distribution  unlimited. 

23 


Set  the  temporary  structures  used  within  this  function. 

/*  For  example,  snd  =  (struct  temporary  *)malloc(sizeof( struct  temporary)); 
and  similarly  for  leveltemp  and  layertemp.  */ 

Initialize  structures  and  temporary  variables  with  the  missing  data/bad  data  indicator.  /*  Missing 
or  bad  data  indicated  by  the  defined  value  ERROR  =  -999.0.  */ 

/*Set  parameters  for  level  and  layer  values.  Here  nht  is  number  of  pressure  levels.  */ 
msg  nht  =  mlevel  nht  -  1 ;  /*  number  of  layer  values  one  less  than  level  values  */ 
size  =  snd  nht; 
msize  =  mlevel  nht; 

pmin  =  log(snd  pressure  at  level[size-l])  -  0.0001; 

/*  Compute  components  from  input  sounding  wind  speed  and  direction.  */ 

start  j  at  -1  /*  Here  j  is  used  to  count  the  number  of  levels  with  wind  data.*/ 

for  i  from  0  to  <  size  incrementing  by  1  /*Here  i  is  the  index  for  the  input  sounding  data.*/ 

{ 

if( input  wind  speed  at  ith  pressure  not  =  ERROR  and  input  wind  direction  at  ith  pressure  not  = 
ERROR) 

{ 

add  1  toj; 

direction  =  -(input  wind  direction  at  pressure  [i])  *  PI/180  +  3*  PI/2; 
snd  u-component  at  pressure[i]  =  cos(direction)  *  input  wind  speed  at  pressure[i]; 
snd  v-component  at  pressure[i]  =  sin(direction)  *  input  wind  speed  at  pressure[i]; 
snd  p[j]  =  log(input  pressure[i]);  /*  Set  to  the  In  of  the  input  pressure.  */ 


} 

wsize  =  j;  /*wsize  is  the  number  of  input  levels  with  wind  data.  */ 

/*  Compute  temperature  and  virtual  temperature.  Temperature  converted  from  C  to  K  as  needed  in 
calling  main  function.  */ 

/*  Compute  level  values.  */ 

/*  Compute  virtual  temperature  from  input  data.  */ 
for  i  from  0  to  <  size  incrementing  by  1 

{  /*  tvfromtemp  computes  virtual  temperature  (Tv)  from  pressure  (P),  sensible  temperature 
(T),  and  relative  humidity  (H)  using  a  standard  method.  Here  pressure  is  the  input  pressure,  not  In 
pressure.*/ 

virtual  temperature  at  input  data  line[i]  of  input  =  tvfromtemp(temperature,  pressure  ,  relative 
humidity)  at  data  line[i]); 

} 


/*  Set  up  temporary  variables  for  use  in  level  and  layer  functions  as  needed.  Height  is  denoted  as 
Z.  */ 

start  j  at  -1  /*  Here  j  is  used  to  count  the  number  of  levels  with  T,  p,  and  H  data.  */ 

for  i  from  0  to  <  size  incrementing  by  1  /*Here  i  is  the  index  for  the  input  sounding  data.*/ 

{ 

if(height,  T,  and  p  at  input  data  level[i]  not  =  ERROR) 

{ 

add  1  toj; 

snd  T[i]  =  input  sounding  T  at  input  data  level[i]; 
snd  Tv[i]  =  input  sounding  Tv  at  input  data  level[i]; 
snd  Z[i]  =  input  sounding  Z  at  input  data  level[i]; 
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snd  H[i]  =  input  sounding  H  at  input  data  level[i]; 
/*  Ln  of  pressure  values  for  snd  computed  earlier.*/ 


} 

} 

tsize  =  j;  /*tsize  is  the  number  of  input  levels  with  T,  Z,  and  H  data.  */ 

for  i  from  0  to  <  msize  incrementing  by  1  /*  Initialize  temporary  level  and  layer  values  of  the  ln 
of  pressure.*/ 

{ 

leveltemp  p[i]  =  upper  boundary  level[i]  p; 
layertemp  p[i]  =  leveltemp  p[i]; 

} 


/*Compute  level  values  for  T,  Tv,  H,  and  Z  where  the  level  function  for  pressure  is  described  in 
this  appendix.*/ 

/*ln  the  function  description  the  respective  variables  or  arrays  are  pmin,  msize,  leveltemp  p,  input 
value,  leveltemp  value,  and  input  p.  */ 

levelprs(pmin,  msize,  leveltemp  p  array,  snd  T  array,  leveltemp  T  array,  snd  p  array); 
levelprs(pmin,  msize,  leveltemp  p  array,  snd  Tv  array,  leveltemp  Tv  array,  snd  parray); 
levelprs(pmin,  msize,  leveltemp  p  array,  snd  H  array,  leveltemp  H  array,  snd  p  array); 
levelprs(pmin,  msize,  leveltemp  p  array,  snd  Z  array,  leveltemp  Z  array,  snd  p  array); 

/*Compute  layer  values  for  T,  Tv,  H,  and  Z  where  the  layer  function  for  pressure  is  described  in 
this  appendix.*/ 

/*ln  the  function  description  the  respective  variables  or  arrays  are  pmin,  size,  htsize,  zh,  value, 
lev_value,  lay_value,  z..  */ 

layerprs(pmin,  tsize,  msize,  layertemp  p  array,  snd  T  array,  leveltemp  T  array,  layertemp  T 
array,  snd  p  array); 

layerprs(pmin,  tsize,  msize,  layertemp  p  array,  snd  Tv  array,  leveltemp  Tv  array,  layertemp  Tv 
array,  snd  p  array); 

layerprs(pmin,  tsize,  msize,  layertemp  p  array,  snd  H  array,  leveltemp  H  array,  layertemp  H 
array,  snd  p  array); 

layerprs(pmin,  tsize,  msize,  layertemp  p  array,  snd  Z  array,  leveltemp  Z  array,  layertemp  Z 
array,  snd  p  array); 

/*  Compute  Wind  Speed  and  Direction  */ 

/*  Compute  level  values  of  wind  components  (u,  v).  In  this  version  p  is  the  ln  of  pressure.  */ 

levelprs(pmin,  msize,  leveltemp  p  array,  snd  u  array,  leveltemp  u  array,  snd  p  array); 
levelprs(pmin,  msize,  leveltemp  p  array,  snd  v  array  ,  leveltemp  v  array,  snd  p  array); 

/*  Compute  level  wind  speed  &  direction  from  components  (u,  v).  */ 

for  i  from  0  to  msize  incrementing  by  1 

{ 

level[i]  value  of  wind  direction  =  (2*  PI  -  atan2(leveltemp  u[i],  -leveltemp  v[i])*  1 80/  PI; 
if(level[i]  direction  >  360) 
subtract  360  from  level[i]  direction; 

level[i]  value  of  wind  speed  =  sqrt(leveltemp  u[i]  *  leveltemp  u[i]  +  leveltemp  v[i]  * 
leveltemp  v[i]); 

} 
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/*  Compute  layer  values  of  components  (u,  v).  */ 

layerprs(pmin,  wsize,  msize,  layertemp  p  array,  snd  u  array,  leveltemp  u  array,  layertemp  u 
array,  snd  p  array); 

layerprs(pmin,  wsize,  msize,  layertemp  p  array,  snd  v  array,  leveltemp  v  array,  layertemp  v 
array,  snd  p  array); 

/*  Compute  message  layer  wind  speed  &  direction  from  components  (u,  v).  Message  values  are 
the  layer  values  in  this  function.  */ 

surface  message  wind  speed  =  input  surface  wind  speed;  /*  i  =0  at  surface*/ 
surface  message  wind  direction  =  input  surface  wind  direction; 
surface  message  u  =  computed  surface  u; 
surface  message  v  =  computed  surface  v; 

for  i  from  1  to  <  msize  incrementing  by  1 

{  /*  Here  message  level[i]  value  is  the  resultant  layer[i]  value.  */ 

message  level[i]  u  =  layertemp  u[i-l]; 
message  level[i]  v  =  layertemp  v[i-l]; 

message  level  [i]  wind  direction  =  (2*  PI  -  atan2(message  level[i]u,  -message  level[i] 
v))*180/PI; 

if  (message  level[i]direction  >  360) 
subtract  360  from  message  levelfi]  direction; 

message  level  [i]  wind  speed  =  sqrt(message  level[i]  u  *  message  level  [i]  u  +  message 
levelfi]  v  *  message  levelfi]  v); 

} 


/*  Load  pressure  values  into  the  message  structures.  */ 

for  i  from  0  to  <  msize  incrementing  by  1 

{ 

message  level[i]  p  =  level[i]  p;  /*  “message  level”  refers  to  resultant  layer.  */ 

I 


/*  Load  T,  Tv,  H,  and  Z  values  into  level  and  message  structures.  */ 

for  i  =  0  to  <  msize  incrementing  by  1  /*  level  values  */ 

{ 

levelfi]  T  =  leveltemp  T[i]; 
level[i]  Tv  =  leveltemp  Tv[i]; 
levelfi]  H  =  leveltemp  H[i]; 
levelfi]  Z  =  leveltemp  Z[i]; 


message  T  at  surface  =  input  surface  T;  /*  Message  values  for  surface,  where  i  =0.  */ 
message  Tv  at  surface  =  Tv  from  surface  data; 
message  H  at  surface  =  input  surface  H; 
message  Z  at  surface  =  input  surface  Z; 

for  i  from  1  to  i  <  msize  incrementing  by  1 

{ 

message  levelfi]  T  =  layertemp  T[i-1];  /*  Message  values  above  surface.  */ 
message  levelfi]  Tv  =  layertemp  Tv[i-1]; 
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message  level [i]  H  =  layertemp  H[i-1]; 
message  level [i]  Z  =  layertemp  Z[i-1]; 

} 


/*  Load  in  site  information  (date,  time,  lat.  Ion,  etc.).  */ 

message  site  information  =  input  site  information; 

/*Free  up  temporary  arrays  in  order  to  release  memory.  */ 

free(snd); 

free(leveltemp); 

free(layertemp); 

return  to  calling  function; 

} 


A-2.2  Alternate  Version  of  msgvaluesprs  Function 

Section  A-2.1  describes  the  primary  version  of  msgvaluesprs.  An  alternate  version 
computes  height  (zcomp)  using  the  hypsometric  fonnula.  The  hypsometric  fonnula 
may  be  found  in  textbooks  and  the  Glossary  of  Meteorology.2  The  main  difference 
from  the  primary  version  of  msgvaluesprs  is  the  use  of  the  procedure  that  uses  the 
hypsometric  fonnula  instead  of  linear  interpolation  in  terms  of  ln(P).  Only  the 
additional  procedure  is  shown  since  the  rest  of  the  function  is  essentially  the  same 
except  for  the  removal  of  statements  for  computation  of  heights  as  a  function  of 
ln(P).  Here  p  is  used  to  represent  ln(P).  As  in  the  previous  section  ‘C’  type  syntax 
is  used  in  the  pseudo  code.  The  natural  logarithm  in  ‘C’  is  written  as  log  vs.  In.  The 
structure  and  variable  names  are  the  same  as  in  the  primary  version  in  the  preceding 
section.  Note  that  snd(p)  was  set  to  snd(ln  of  pressure)  prior  to  this  excerpt  from 
the  alternate  version  of  msgvaluesprs. 

/*  Start  of  height  computation  section.  P  is  pressure  and  p  is  ln(P).*/ 
surface  level  Z  =  snd  surface  Z; 


j=o; 

for(i=l;  i<msize;  i++)  /*Values  for  temporary  sound  structure  set  before  level 

calculations  above,  “mlevel”  refers  to  pressure  level  values.  */ 

{ 

while(snd  p[j]  >  log( mlevel  P  at  level[i]) 

{ 

j++; 

} 


if(snd  p  [j]  =  log( mlevel  P  at  level[i]))  /*Sounding  P  =  upper  boundary  P*/ 


2 

“  Glickman  TS,  managing  editor.  Glossary  of  meteorology.  2nd  ed.  Boston  (MA):  American 
Meteorological  Society;  2000. 
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{ 

mlevel  Z  at  level[i]  =  snd  Z[j] ; 

} 

else 

{ 

if(snd  p[j-l]  >  log(mlevel  P  at  level[i-l])  &&  j  >  0)  /* Sounding  P  higher  or  height  of 
pressure  level  lower.*/ 

{ 

mlevel  Z  at  level[i]  =  zcomp(mlevel  Tvat  level[i].,  mlevel  Tv  at  level[i-l], 

mlevel  P  at  level[i-l],  mlevel  P  at  level[i],  mlevel  Z  at  level[i-l]); 

} 

else 

{ 

mlevel  Z  at  level[i]  =  zcomp(mlevel  Tvat  level[i],  snd  Tv[j-1], 

exp(snd  p[j- 1]),  mlevel  P  at  level[i],  snd  Z[j- 1]);  /*  snd  p  is  ln(input  P)*/ 

} 

} 

} 

/*  End  of  computation  of  level  pressures.*/ 

/*  Compute  layer  height  values,  msg  values  =  values  for  surface  +  values  for  layers.*/ 
msg  Z  at  surface  =  mlevel  Z  at  surface;  /*  surface  =  level[0]  */ 

j=i; 

for(i=l  ;i<msize;i++) 

{ 

while(snd  p[j]  >  (log(mlevel  P  at  level[i])  +  log(mlevel  P  at  level[i-l]))*0.5) 

{ 

j++; 

} 

if(snd  p[j]  ==  (log( mlevel  P  at  level[i])  +  log(mlevel  P  at  level[i-l]))*0.5) 

{ 

msg  Z  at  level[i]  =  snd  Z[j]; 

} 

else 

{ 

if(snd  p[j-l]  <  (log(mlevel  P  at  level[i-l])  +  log( mlevel  P  at  level[i-2]))*0.5) 

! 

Tv  =  msg  Tv  at  level[i]; 

Tvo  =  mlevel  Tv  at  level[i-l]; 

P  =  (mlevel  P  at  level[i]  +  mlevel  P  at  level[i-l])*0.5; 

Po  =  mlevel  P  at  level[i-l]; 

Zo  =  mlevel  Z  at  level[i-l]; 

msg->level[i].hgt  =  zcomp(Tv,  Tvo,  Po,  P,  Z0); 

} 

else 

{ 
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Tv  =  msg  Tv  at  level[i]; 

Tvo  =  snd  Tv[j]; 

P  =  (mlevel  P  at  level[i]  +  mlevel  P  at  level[i-l])*0.5; 
Po  =  exp(sndp[j-l]); 

Z0  =  snd  Z[j-1]; 

msg  Z  at  level[i]  =  zcomp(Tv,  Tvo,  Po,  P,  Zo); 


/*  End  of  computation  of  layer  pressure  values.*/ 


A-2.3  Interpolation  and  Integrated  Mean  Functions 

The  levelprs  and  layerprs  functions  are  close  to  the  level  and  layer  functions 
described  in  Cogan1  that  perform  linear  interpolation  and  computation  of  integrated 
mean  layer  values  based  on  height.  The  main  differences  are  the  substitution  of  a 
pressure  variable  (p)  for  height,  interpolation  and  computation  of  integrated  mean 
from  highest  to  lowest  p  values  vs.  lowest  to  highest  height,  and  a  procedure  to 
check  that  p  is  greater  or  equal  to  some  minimum  value,  pmin,  vs.  less  than  or  equal 
to  a  maximum  height,  zmax.  The  variable  p  in  these  functions  may  represent  P  or 
ln(P),  and  similarly  for  pi.  As  in  the  previous  section  ‘C’  type  syntax  is  used  in  the 
pseudo  code. 

A-2.3.1  Description  of  levelprs 

Pass  in  the  minimum  value  of  p  (pmin),  the  number  of  the  user-defined  pressure  levels  (prsize), 
the  user-defined  pressure  levels  (pi  array),  the  input  values  of  a  variable  (value  array),  and  the 
input  p  data  (p  array).  Pass  out  of  the  function  the  level  values  of  the  variable  (lev  value  array). 

{ 

set  i=0  and  j=0;  /*  Starting  at  surface  or  level  nearest  surface.*/ 
while(p[j]  >  pmin  and  pl[i]  >  pmin  and  i  <  prsize) 

{ 

if(pl[i]  less  than  or  equal  to  p[0])  /*  p[0]  is  the  surface  value  (or  for  p  level  nearest  the  surface) 

*/ 

{ 

if(pl[i]  equals  p[j])  /*  Input  pressure  =  defined  pressure  level.*/ 

{ 

lev_value[i]  =value[j]; 
add  1  toj; 

} 

else 

{ 

if(p[j]  less  than  pl[i])  /*  use  to  go  to  previous  input  p  (except  for  surface  or  level  nearest 
surface).*/ 

{ 

subtract  1  fromj; 

} 

lev  valuefi]  =  valuefj]  -  (valuefj]  -  value[j+l])*(p[j]-pl[i])/(p[j]-p[j+l]); 
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/interpolate  to  obtain  levvalue  (value  for  given  pressure  level).*/ 


} 

} 

while(pl[i+l]  less  than  or  equal  to  p[j+l]  and  p[j+l]  greater  than  or  equal  to  pmin  and  p[j]  >  - 
999)  /*-999  indicates  no  data  value*/ 

{ 

add  1  toj; 

} 

add  1  to  i;  /*  Go  to  next  input  p  level.  */ 


return; 

}  /*  End  of  level  function.  */ 


A-2.3.2  Description  of  layerprs 

Pass  in  the  minimum  value  of  p  (pmin),  the  number  of  input  data  values  of  the  variable  (size),  the 
number  of  the  user-defined  pressure  levels  (psize),  the  user-defined  pressure  levels  (pi  array),  the 
input  values  of  a  variable  (value  array),  the  level  values  of  the  variable  (lev_value  array),  and  the 
input  p  data  (p  array).  Pass  out  of  the  function  the  layer  values  of  the  variable  (lay_value  array). 

{ 

Set  up  the  temporary  arrays:  tempval,  tempp,  and  mean. 

i  =  0; 

while(i  is  less  than  (psize-1)  and  pl[i+l]  greater  than  or  equal  to  pmin)  /*  Find  layer  means.*/ 

{ 

add  1  to  i; 

set  ind  =  0  and  sum  =  0; 

tempval[ind]  =  lev_value[i-l];  /*  Lower  (higher  p)  boundary  values.*/ 
tempp[ind]  =pl[i-l]; 

for  (j  starting  at  0,  to  the  highest  value  <  size,  incrementing  by  1) 

{ 

if(p[j]  greater  than  pl[i]  and  p[j]  less  than  pl[i-l])  /*  Values  within  layer.*/ 

! 

add  1  to  ind; 
tempval[ind]  =value[j]; 
tempp  [ind]  =  p[j]; 


} 

add  1  to  ind;  /*  Upper  (lower  p)  boundary  level.*/ 

tempval[ind]  =  lev_value[i]; 
tempp[ind]  =  pl[i] ; 

for  (j  starting  at  1,  up  to  and  including  ind,  incrementing  by  1)  /*  Sub-layer  average.*/ 

{ 

mean[j-l]  =  (tempval[j]  +  tempval[j-l])  *  0.5; 

} 


for(j  starting  at  1,  up  to  and  including  ind,  incrementing  by  1)  /*  Proportional  weighting  of 
each  layer.*/ 
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{ 

add  (mean[j-l]  *  (tempp[j-l]  -  tempp[j]))  to  sum; 

} 

lay_value[i-l]  equals  sum  /(pl[i- 1  ]  -  pl[i]);  /*  Mean  layer  value  =  sum/layer  thickness.*/ 

}  /*  end  of  while  loop.*/ 

/*  Free  temporary  arrays.  ******/ 

free(tempval); 

free(tempp); 

free(mean); 

return; 

}  /*  End  of  layer  function.  */ 
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Intentionally  Left  Blank 
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Appendix  B.  Charts  of  Mean  Difference  (MD),  Mean  Absolute 
Difference  (MAD),  Standard  Deviation  (SD),  and  Root  Mean 
Square  Difference  (RMSD)  for  All  Variables 
for  Each  of  the  4  Methods 
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Fig.  B-l  Statistical  values  (MD,  MAD,  SD,  and  RMSD)  for  the  several  variables  for  the  user- 
defined  height  levels  along  with  the  number  of  samples  for  each  height  level 
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Temperature  Difference 
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Fig.  B-l  Statistical  values  (MD,  MAD,  SD,  and  RMSD)  for  the  several  variables  for  the  user- 
defined  height  levels  along  with  the  number  of  samples  for  each  height  level  (continued) 
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Fig.  B-2  Statistical  values  (MD,  MAD,  SD,  and  RMSD)  for  the  several  variables  for  height 
layers  defined  by  the  user  input  height  levels  along  with  the  number  of  samples  for  each  height 
layer 
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Fig.  B-2  Statistical  values  (MD,  MAD,  SD,  and  RMSD)  for  the  several  variables  for  height 
layers  defined  by  the  user  input  height  levels  along  with  the  number  of  samples  for  each  height 
layer  (continued) 
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Fig.  B-3  Statistical  values  (MD,  MAD,  SD,  and  RMSD)  for  the  several  variables  for  pressure 
levels  defined  by  the  user  along  with  the  number  of  samples  for  each  level 
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Fig.  B-3  Statistical  values  (MD,  MAD,  SD,  and  RMSD)  for  the  several  variables  for  pressure 
levels  defined  by  the  user  along  with  the  number  of  samples  for  each  level  (continued) 
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List  of  Symbols,  Abbreviations,  and  Acronyms 


AGL 

above  ground  level 

ECMWF 

European  Center  for  Medium-Range  Weather  Forecasting 

GFS 

Global  Forecast  System 

IMS 

Ice  Mapping  System 

JMA 

Japan  Meteorological  Agency 

M 

mean  value 

MAD 

mean  absolute  difference 

MADIS 

Meteorological  Assimilation  Data  Ingest  System 

MAE 

mean  absolute  error 

METAR 

Meteorological  Aerodrome  Report 

MD 

mean  difference 

MM  5 

Meteorological  Model  Fifth  Generation 

MSL 

mean  sea  level 

MYJ 

Mellor-Yamada-Janjic  scheme 

NCAR 

National  Center  for  Atmospheric  Research 

NCEP 

National  Center  for  Environmental  Prediction 

NCL 

National  Command  Fanguage 

NOHRSC 

National  Operational  Hydrologic  Remote  Sensing  Center 

NWP 

numerical  weather  prediction 

PWV 

Precipitable  Water  Vapor 

RAOB 

radiosonde  observation 

RASS 

radio  acoustic  sounding  systems 

RMSD 

root  mean  square  difference 

RMSE 

root  mean  square  error 

SAO 

Surface  Aviation  Observation 

SD 

standard  deviation 
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SNODAS  Snow  Data  Assimilation  System 


WMO  World  Meteorological  Organization 

WRF  Weather  Research  and  Forecasting 
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